package com.fy.utils;

import org.apache.commons.math3.distribution.TDistribution;
import org.apache.commons.math3.stat.descriptive.moment.StandardDeviation;

/**
 * 统计学 p值 T值 方差计算
 */
public class MatlabUtil {
    private double[] x;
    private double[] y;
    StandardDeviation standardDeviation =new StandardDeviation();


    public MatlabUtil() {

    }

    public MatlabUtil(double[] x, double[] y) {
        this.x = x;
        this.y = y;
    }

    public int getXSize() {
        return x==null?0:x.length;
    }
    public int getYSize() {
        return y==null?0:y.length;
    }

    /**
     * 平均值
     * @return
     */
    public double getXMean() {
        int n=x.length;
        double sum=0;
        for (double d : x) {
            sum=sum+d;
        }
        return sum/n;
    }

    /**
     * 平均值
     * @return
     */
    public double getYMean() {
        int n=y.length;
        double sum=0;
        for (double d : y) {
            sum=sum+d;
        }
        return sum/n;
    }

    /**
     * 标准差
     * @param x
     * @return
     */
    public double getStandard(double[] x) {
        return standardDeviation.evaluate(x);
    }

    /**
     * T值
     * @return
     */
    public double calculateTvalue() {
        double a=getXMean()-getYMean();
        double q1=getStandard(x);
        double q2=getStandard(y);
        double s1=q1*q1;
        double s2=q2*q2;
        double b=Math.sqrt(s1/x.length+s2/y.length);
        return a/b;
    }

    /**
     * 自由度 N
     * @return
     */
    public double getDegreesOfFreedom() {
        double q1=getStandard(x);
        double q2=getStandard(y);
        double s1=q1*q1;
        double s2=q2*q2;
        double a=(s1/x.length+s2/y.length)*(s1/x.length+s2/y.length);
        double b=((s1/x.length)*(s1/x.length))/(x.length-1)+((s2/y.length)*(s2/y.length))/(y.length-1);
        Double free=a/b;
        return free.doubleValue();
    }

    /**
     * T值
     * @return
     */
    public double getPValue() {
        double free=getDegreesOfFreedom();
        double t=calculateTvalue();
        TDistribution td=new TDistribution(free);
        double cumulative = td.cumulativeProbability(t);
        double p;
        if(t>0) {
            p=(1-cumulative)*2;
        }else {
            p=cumulative*2;
        }
        return p;
    }

    /**
     *T值
     * @param xmean 样本1均值
     * @param ymean 样本2均值
     * @param xq 样本1标准差
     * @param yq 样本2标准差
     * @param xcount 样本1数量
     * @param ycount 样本2数量
     * @return
     */
    public double calculateTvalue(double xmean,double ymean, double xq,double yq,double xcount,double ycount) {
        double a = xmean-ymean;
        double q1=xq;
        double q2=yq;
        double s1=q1*q1;
        double s2=q2*q2;
        double b=Math.sqrt(s1/xcount+s2/ycount);
        return a/b;
    }

    /**
     * 自由度 N值
     * @param xq
     * @param yq
     * @param xcount
     * @param ycount
     * @return
     */
    public double getDegreesOfFreedom(double xq,double yq,double xcount,double ycount) {
        double q1=xq;
        double q2=yq;
        double s1=q1*q1;
        double s2=q2*q2;
        double a=(s1/xcount+s2/ycount)*(s1/xcount+s2/ycount);
        double b=((s1/xcount)*(s1/xcount))/(xcount-1)+((s2/ycount)*(s2/ycount))/(ycount-1);
        Double free=a/b;
        return free.doubleValue();
    }

    /**
     * P值
     * @param xmean
     * @param ymean
     * @param xq
     * @param yq
     * @param xcount
     * @param ycount
     * @return
     */
    public double getPValue(double xmean,double ymean, double xq,double yq,double xcount,double ycount) {
        double free= getDegreesOfFreedom(xq,yq,xcount,ycount);
        double t = calculateTvalue(xmean,ymean,xq,yq,xcount,ycount);
        TDistribution td=new TDistribution(free);
        double cumulative = td.cumulativeProbability(t);
        double p;
        if(t>0) {
            p=(1-cumulative)*2;
        }else {
            p=cumulative*2;
        }
        return p;
    }

    public static void main(String[] args) {
        MatlabUtil matlabUtil = new MatlabUtil();
        System.out.println("T:"+ matlabUtil.calculateTvalue(52,50,4,4,10,8));
        System.out.println("N:"+ matlabUtil.getDegreesOfFreedom(4,4,10,8));
        System.out.println("P:"+ matlabUtil.getPValue(52,50,4,4,10,8));
    }
}
